Phase transition in a spatial Lotka-Volterra model 
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Spatial evolution is investigated in a simulated system of nine competing and mutating bacterium strains, 
which mimics the biochemical war among bacteria capable of producing two different bacteriocins (toxins) at 
most. Random sequential dynamics on a square lattice is governed by very symmetrical transition rules for 
neighborhood invasions of sensitive strains by killers, killers by resistants, and resistants by sensitives. The 
community of the nine possible toxicity /resistance types undergoes a critical phase transition as the uniform 
transmutation rates between the types decreases below a critical value Pc above which all the nine types of 
strain coexist with equal frequencies. Passing the critical mutation rate from above, the system collapses into 
one of three topologically identical (degenerated) states, each consisting of three strain types. Of the three 
possible final states each accrues with equal probability and all three maintain themselves in a self-organizing 
polydomain structure via cyclic invasions. Our Monte Carlo simulations support that this symmetry-breaking 
transition belongs to the universality class of the three-state Potts model. 

PACS numbers: 87.23. Cc, 05.10.-a, 05.40. Fb, 64.60.Ht 



Many species of bacteria have recently been shown to 
excrete toxic subtances that are very effective against 
strains of the same or closely related species not pro- 
ducing the corresponding resistance factor With 
respect to a certain toxin a species consists of colonies 
of three possible types: Sensitive (5), Killer (K), or Re- 
sistant (i?). Killer strains produce the toxin and a resis- 
tance factor that prevents suicide; resistant strains only 
produce the resistance factor, and sensitives produce nei- 
ther. An S colony can always be invaded and displaced 
hy a, K propagulum because K can kill S using the toxin. 
A K colony can be invaded and ultimately displaced by 
an R propagulum because the resistant type is immune to 
the toxic effect and it does not carry the metabolic bur- 
den of synthetizing the toxin, thus it achieves a higher 
growth rate and competitive dominance over K. The 
sensitive (S) type can in turn displace the resistant by 
competition, for it does not even pay the metabolic cost 
of producing the resistance factor. The resulting cyclic 
pattern of competitive dominance {K beats S beats R 
beats K) is a striking realization of the well-known Rock- 
Scissors-Paper game by a biological entity. Other 
cyclic dominance systems are almost unknown in ecol- 
ogy, but given the extraordinary significance of bacterial 
communities in virtually all ecological systems, the prob- 
lem is well worth detailed theoretical studies. 

Some theoretical aspects of cyclic dominance have al- 
ready been thoroughly investigated . In the simplest 
spatial (lattice) version of a cyclic Lotka-Volterra system 
|,| the species are distributed on a d-dimensional lat- 
tice, and invasions are confined to nearest neighbor sites 
with uniform rates. |^,^ Analytical and numerical cal- 
culations have proven that fixation occurs if the number 



of species exceeds a critical value dependent on dimen- 
sion d, otherwise a self-organizing domain structure is 
maintained for d > 2. which comprises rotating vor- 
tices and antivortices in three-species models. The 
present work is meant to demonstrate that extending the 
cyclic dominance approach to a two-toxin bacterial com- 
munity with mutation results in a remarkable enrichment 
of interesting dynamical phenomena, compared to what 
is already known. 

The relevant biological details of bacteriocin systems 
are the following: The genes coding for the toxin and 
the resistance factor are usually both sitting on an extra- 
chromosomal DNA-ring in the citoplasm (called plasmid) 
that the bacterium can lose and obtain without any im- 
mediate deleterious impact. Each of the two genes on the 
plasmid can be switched off by DNA mutation. Thus all 
possible mutational transformations are possible in prin- 
ciple, but — supposing the mutant does not disperse far 
immediately — those having a visible effect are only the 
ones after which the mutant defeats the resident strain 
from which it emerged. Obviously, S — > K, K R, 
and R ^ S are mutations followed by competitive dis- 
placement of the resident, i.e., they are permitted, but 
the reverse ones are immediately eliminated by the resi- 
dent population. S ^ K involves obtaining a complete 
plasmid which is possible, e.g., through genetic transfor- 
mation or a sexual event called conjugation; the other 
two viable mutations are realized by switching off the 
toxin gene and the resistance gene, respectively, on an 
existing plasmid. 

Most bacteria are capable of producing more than one 
toxin and/or the corresponding resistance factors simul- 
taneously. If the maximum number of toxin types is two. 
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the number of possible toxicity/resistance combinations 
in a strain is nine. These are: SS, SK, SR, KS, KK, 
KR, RS, RK, and RR. Here we confine our attention 
to this two-toxin system, denoting the actual state of a 
bacterium colony by an index number from to 8 in the 
order above. The topology of the dominance relations 
among the states is illustrated in Fig. |l|. The biological 
justification for this topology is straightforward: double 
dominance of strain A above strain B means that A har- 
bors dominant genes on both plasmids compared to B; 
single dominance means that one gene of A is dominant, 
the other is identical to that of the correspondig gene in 
B; no dominance follows either if the corresponding genes 
are both identical, or if the two genes play draw (i.e., A 
wins with one gene, and B with the other). 




FIG. 1. Topology of dominance in the nine-species model. 
The single and double lines with arrow indicate single and 
double dominances as described in the text. 

We use a square lattice of size L x L with periodic 
boundary conditions as the arena for interaction. Each 
lattice site i is occupied by a colony of one of the nine tox- 
icity/resistance types (we call them "species" for brevity 
in the sequel) . Assignment of a certain value of the state 
variable = 0, . . . , 8 to a certain site indicates the pres- 
ence of the corresponding species on that site. The simu- 
lation works by iterating the following elementary steps: 
At a randomly chosen site one of the two possible mutants 
replaces the resident colony with probability P, otherwise 
the resident colony fights with a randomly chosen nearest 
neighbor colony by mutually invading propagules (with a 
probability 1 — 2P). The outcome of the battle between 
two neighbors (i and j) depends on the dominance rela- 
tions between them: i displaces j and takes over its site 
if Si is dominant over Sj (cf. Fig. ^ . If the neighbors are 
equivalent or neutral to each other (play draw), nothing 
happens. 

Lattice size varies from L — 400 to 3000 in the simu- 
lations. At t = the species were distributed at random 



with uniform probabilities on the lattice in all simulation 
runs. The control parameter was mutation rate, vary- 
ing from to 1/2. We have recorded the time series of 
species concentrations and correlation functions, the lat- 
ter of which served as the basis for calculating correlation 
length. After a suitable thermalization time we have av- 
eraged these data over some sampling time chosen to be 
long enough for providing sufhcient accuracy in spite of 
the occassionally very high concentration fluctuations. 

In no-mutation runs (P = 0) we have observed an in- 
teresting domain size increase phenomenon in the time 
series of spatial patterns which develop. One can dis- 
tinguish three equivalent types of growing domains con- 
sisting of the 0+4+8, 1+5+6, and 2+3+7 species re- 
spectively. Inside these three domains a self-organizing 
structure is maintained through the mechanism described 
by Tainaka for the simplest 3-species cyclic dominance 
model 1^. We call these domains "alliances" henceforth, 
but note that the species within an alliance are in fact 
the worst enemies: each alliance consists of species cycli- 
cally double-dominating each other. Our reason for this 
terminology will be clear in a minute. 

The three alliances are given approximately equal ter- 
ritories at start, but the system slowly drifts towards a 
single-domain state in all simulations. Each alliance has 
the same chance to take over. 

Alliances defend themselves against the external inva- 
sion of "alien" species with a peculiar mechanism. One 
can easily check that any external invador can attack 
only one of the species within an alliance, and the in- 
vador is eliminated from the domain of the alliance by 
the species actually controlling the attacked one within 
the alliance. This means that the invador is wiped away 
by the toughest within-domain enemy of the attacked 
species, thus maintaining the self-organizing structure 
and the integrity of the domain with the very same mech- 
anism. 

The self-protection of alliances against external in- 
vadors can also be observed for small mutation rates as 
illustrated in Fig. |^. In this snapshot similar symbols are 
used for species belonging to the same alliance. Namely, 
different strip widths (or box sizes) distinct the species 
within the three alliances represented by horizontal and 
vertical strips and closed squares respectively. This figure 
illustrates that the mutants and their offspring can form 
only small, temporary islands in the sea of the dominant 
domain (0+4+8) represented by vertical strips. Clearly, 
the concentrations of minority species increase with mu- 
tation rate P. 

The average concentrations of the species become equal 
if the mutation rate exceeds a critical value Pc- This con- 
tinuous transition is accompanied by a divergence in both 
the fluctuations and the correlation length. A similar 
phase transition occurs in the three-state Potts model 
||ll|,|2|, therefore we have adopted the numerical tech- 
niques suggested by Binder ||l^ for the quantitative anal- 
yses. 
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FIG. 2. Snapshot of species pattern for P = 0.0003. The 
different symbols magnified at the top represent the species 
s = 0, . . . , 8 from left to right. 

In order to reduce relaxation time in systematic inves- 
tigations, the random initial state contained only species 
0, 4, and 8 at small values of P. The MC simulations were 
performed on large lattices (L > 1500) and long sampling 
times (t > 3 • 10^ MC steps per sites) in the vicinity of 
the critical point. For such large lattice sizes the domi- 
nance of the 0-1-4+8 alliance was maintained during the 
simulations at all tested values of P below R,- 
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FIG. 3. Average concentrations of species as a function of 
mutation rate. The inset shows the log-log plot of the order 
parameter vs. Pc — P- 



dicated in Fig. ^. Statistical errors are small, comparable 
to the line thickness of the figure, except for the very close 
vicinity of the critical point. 

Figure ^ shows that the P-dependence of average con- 
centrations can be characterized by a single order param- 
eter m as 
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Co = C4 = C8 = - -I- m , 



(1) 
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Ci = C2 = C3 = C5 = C6 = C7 = - - -m . 



According to our MC simulations, m follows a power law 
behavior in the close vicinity of Pc, i.e. 



m oc (Pc - P) 
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(2) 



if P < Pc (see the inset in Fig. whereas the order 
parameter remains zero for P > Pc- Numerical fitting 
yields (3 = 0.110(5) and P^ = 0.0004333(5). This value of 
the (3 exponent is in good agreement with the theoretical 
prediction (/3 — 1/9) obtained for the three-state {Q = 3) 
Potts model |12|. This is not surprising, given that a 
large class of two-state dynamical systems exhibits phase 
transition belonging to the universality class of the Ising 
model [|l^ and the Q-state Potts model was introduced 
as a generalization of the Ising model ||ll|,|l^ . 

To obtain further evidence, we have also studied some 
other quantities characterizing the critical behavior. For 
example, the fluctuation of the order parameter defined 
as X = N{{m — (m))^) can be well approximated by a 
power law; X ~ 1^ ^ ^c|^ in the vicinity of Pc- Be- 
low and above the critical point, numerical fitting yields 
7 = 1.3(2) and 7' = 1.43(4) respectively, which values 
agree with the theoretical prediction (7 = 7' = 13/9) 
p2| . The investigation of the cumulant of the order pa- 
rameter |l^ for small lattice sizes (L = 60, 100, 200) 
supports the presence of a continuous transition at the 
critical mutation rate Pc- Furthermore, we have deter- 
mined the correlation function C{x) for P > Pc, which 
characterizes the probability of finding the same species 
on two sites at a distance x away from each other. In the 
vicinity of Pc, two different characteristic lengths can be 
obtained from C{x) (see the inset in Fig. The short- 
est correlation length is proportional to the linear size of 
a domain within the alliance, and this quantity remains 
finite if P —> Pc. The longest correlation length is more 
interesting, because it characterizes the linear size of the 
alliance and diverges if P — > Pc- More precisely, this 
correlation length can be well described by a power law 
(C ^ (P^PcT) as illustrated in Fig. ^. Numerical fitting 
predicts 1/ = 0.82(4), in agreement with the theoretical 
prediction v ^ 5/6 for the two-dimensional, three-state 
Potts model [HI. 



The simulations suggest a continuous transition as in- 
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FIG. 4. Log-log plot of the correlation length characteristic 
to the linear extension of alliances as a function of P — Pc- 
The solid line represents the fitted power law behavior with 
a slope of -0.82. The arrow points to the value of the correla- 
tion length characteristic to the linear domain size inside an 
alliance for P = 0. The inset shows the lin-log plot of the 
correlation function for P — 0.00055. 

We have also investigated the model analytically, by 
evaluating the probability of configurations on two adja- 
cent sites (pair approximation method). Even though the 
number of possible pair configurations is large (9^), con- 
sidering the (translation, rotation, reflection, and cyclic) 
symmetries of the system the number of different pair 
probabihties reduces to as few as four. For sufficiently 
large P, this method gives a good approximation for the 
behavior of the simulation model (i.e., it predicts van- 
ishing correlations if P —>■ 1/2). However, it shows no 
sign of the phase transition found in the MC simulations 
at Pc- This failure of the pair approximation method in 
showing the phase transition is related to the key role 
that interfacial invasion plays in the development of the 
self-organizing domain structure. [[7|p0|] This feature lim- 
its the techniques we can use for further investigations. 

In conclusion, our MC simulations justify that the 
present model exhibits a critical phase transition accom- 
panied with spontaneous symmetry-breaking, in close 
analogy to the well-known Potts model. Here the muta- 
tion rate P plays the role of the control parameter whose 
increase drives the system towards the symmetric sta- 
tionary state in which all the nine species are present 
with the same probability. Conversely, for low muta- 
tion rate P the system drifts towards the dominance of 
an alliance consisting of three species whose survival is 
maintained by cyclic within-domain invasion. Due to the 
very symmetric topology of species dominance relations 
the system admits three equivalent alliances. Numeri- 
cal analysis of the critical behavior (/3, v, and 7 expo- 
nents) supports our conjecture that the observed phase 
transition belongs to the universality class of the three- 



state Potts model. The most surprising result of this 
work is that cyclic invasion is capable of providing pro- 
tection (stability) for alliances consisting of mortal ene- 
mies (species double-dominating each other) under some 
particular conditions hidden in the topology of the inter- 
action. Further systematic research is required to clar- 
ify the conditions for the emergence of such defensive 
alliances accompanied by a reduction in the number of 
surviving species for more general interaction topologies. 
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